PRO FITS
; Read FITS files. By Jianfeng Zhou
; Add more comments.

C=[0,1,-1,-1,1,-1,1,1,-1,1,1,1,-1,1,-1,$
   1,1,1,-1,-1,-1,-1,-1,-1,1,1,-1,-1,1,$
   1,-1,-1,-1,-1,-1,-1,1,1,1,-1,1,-1,1,$
   1,1,-1,1,1,-1,1,-1,-1,1]

A=INTARR(53,53)
A=(C#C+1)/2
A[*,0]=0
A[0,1:52]=1

M=INTARR(95,95)
M[0:47,0:47]=ROTATE(A[0:47,0:47],2)
M[48:94,0:47]=ROTATE(A[1:47,0:47],7)
M[0:47,48:94]=ROTATE(A[0:47,1:47],5)
M[48:94,48:94]=A[1:47,1:47]

G=M+M-1
G[47,47]=1

L1=11.2
L2=4.6
G1=INTARR(232,232)
FOR j=116,230 DO BEGIN
 FOR i=116,230 DO BEGIN
  X1=FLOOR((L2*(i-116)+L1/2)/L1)
  X2=FLOOR((L2*(i-116)+L1/2+1)/L1)
  Y1=FLOOR((L2*(j-116)+L1/2)/L1)
  Y2=FLOOR((L2*(j-116)+L1/2+1)/L1)
  CASE 1 of
   (X1 EQ X2) && (Y1 EQ Y2): BEGIN
     G1[j,i]=G[47+Y1,47+X1]
     END
   (X1 EQ X2) && (Y1 NE Y2): BEGIN
     G1[j,i]=G[47+Y1,47+X1]*(Y2*L1-L2*(j-116))/L2+G[47+Y2,47+X1]*(L2*(j-116+1)-Y2*L1)/L2
     END
   (X1 NE X2) && (Y1 EQ Y2): BEGIN
     G1[j,i]=G[47+Y1,47+X1]*(X2*L1-L2*(i-116))/L2+G[47+Y1,47+X2]*(L2*(i-116+1)-X2*L1)/L2
     END
  ELSE : BEGIN
     G1[j,i]=G[47+Y1,47+X1]*((X2*L1-L2*(i-116))*(Y2*L1-L2*(j-116)))/(L2*L2)+G[47+Y1,47+X2]*((L2*(i-116+1)-X2*L1)*(Y2*L1-L2*(j-116)))/(L2*L2)$
             +G[47+Y2,47+X1]*((X2*L1-L2*(i-116))*(L2*(j-116+1)-Y2*L1))/(L2*L2)+G[47+Y2,47+X2]*((L2*(i-116+1)-X2*L1)*(L2*(j-116+1)-Y2*L1))/(L2*L2)
     END
  ENDCASE
 ENDFOR
ENDFOR  
G1[231,231]=G[94,94]
FOR j=116,230 DO BEGIN
  Y1=FLOOR((L2*(j-116)+L1/2)/L1)
  Y2=FLOOR((L2*(j-116)+L1/2+1)/L1)
  IF Y1 EQ Y2 THEN G1[j,231]=G[47+Y1,94] $
  ELSE G1[j,231]=G[47+Y1,94]*(Y2*L1-L2*(j-116))/L2+G[47+Y2,94]*(L2*(j-116+1)-Y2*L1)/L2
ENDFOR
G1[231,116:230]=G1[116:230,231]
G1[0:115,0:115]=ROTATE(G1[116:231,116:231],2)
G1[116:231,0:115]=ROTATE(G1[116:231,116:231],7)
G1[0:115,116:231]=ROTATE(G1[116:231,116:231],5)

M1=INTARR(232,232)
FOR j=116,230 DO BEGIN
 FOR i=116,230 DO BEGIN
  X1=FLOOR((L2*(i-116)+L1/2)/L1)
  X2=FLOOR((L2*(i-116)+L1/2+1)/L1)
  Y1=FLOOR((L2*(j-116)+L1/2)/L1)
  Y2=FLOOR((L2*(j-116)+L1/2+1)/L1)
  CASE 1 of
   (X1 EQ X2) && (Y1 EQ Y2): BEGIN
     M1[j,i]=M[47+Y1,47+X1]
     END
   (X1 EQ X2) && (Y1 NE Y2): BEGIN
     M1[j,i]=M[47+Y1,47+X1]*(Y2*L1-L2*(j-116))/L2+M[47+Y2,47+X1]*(L2*(j-116+1)-Y2*L1)/L2
     END
   (X1 NE X2) && (Y1 EQ Y2): BEGIN
     M1[j,i]=M[47+Y1,47+X1]*(X2*L1-L2*(i-116))/L2+M[47+Y1,47+X2]*(L2*(i-116+1)-X2*L1)/L2
     END
  ELSE : BEGIN
     M1[j,i]=M[47+Y1,47+X1]*((X2*L1-L2*(i-116))*(Y2*L1-L2*(j-116)))/(L2*L2)+M[47+Y1,47+X2]*((L2*(i-116+1)-X2*L1)*(Y2*L1-L2*(j-116)))/(L2*L2)$
             +M[47+Y2,47+X1]*((X2*L1-L2*(i-116))*(L2*(j-116+1)-Y2*L1))/(L2*L2)+M[47+Y2,47+X2]*((L2*(i-116+1)-X2*L1)*(L2*(j-116+1)-Y2*L1))/(L2*L2)
     END
  ENDCASE
 ENDFOR
ENDFOR  
M1[231,231]=M[94,94]
FOR j=116,230 DO BEGIN
  Y1=FLOOR((L2*(j-116)+L1/2)/L1)
  Y2=FLOOR((L2*(j-116)+L1/2+1)/L1)
  IF Y1 EQ Y2 THEN M1[j,231]=M[47+Y1,94] $
  ELSE M1[j,231]=M[47+Y1,94]*(Y2*L1-L2*(j-116))/L2+M[47+Y2,94]*(L2*(j-116+1)-Y2*L1)/L2
ENDFOR
M1[231,116:230]=M1[116:230,231]
M1[0:115,0:115]=ROTATE(M1[116:231,116:231],2)
M1[116:231,0:115]=ROTATE(M1[116:231,116:231],7)
M1[0:115,116:231]=ROTATE(M1[116:231,116:231],5)

D=READFITS('/workroom2/odyssey/Integral/vela/obs/008600440140/scw/008600440140.001/isgri_cor_shad.fits',EXTEN_NO = 2)
D0=INTARR(361,365)
D0[0:129,0:133]=D

M10=INTARR(463,463)
M10[0:231,0:231]=M1
MG=CONVOL(M10,G1,/EDGE_ZERO,CENTER=0)

S1=CONVOL(D0,G1,CENTER=0,/EDGE_ZERO)

CLEAN=INTARR(361,365)
CC=0L
REPEAT BEGIN
P=S1[0,0]
U=0
V=0
DIF=INTARR(361,365)
FOR i=0,364 DO BEGIN
 FOR j=0,360 DO BEGIN
  IF S1[j,i] GT P THEN BEGIN
    P=S1[j,i]
    U=i
    V=j
  ENDIF 
 ENDFOR
ENDFOR
DIF[V,U]=0.03*P
S10=INTARR(823,827)
S10[231:591,231:595]=S1
MG0=INTARR(823,827)
MG0[V:V+462,U:U+462]=0.03*MG
S10=S10-MG0
S1=S10[231:591,231:595]
CLEAN=CLEAN+DIF
CC=CC+1
;R=MAX(IMSL_COVARIANCES(S1))
ENDREP UNTIL CC GT 20000;(R LT 1000000) || (CC GT 20000)
PRINT,CC;,R
S1=CLEAN+S1


WINDOW,0,Xsize=361,Ysize=365
TVSCL,S1
ZOOM,FACT=2,XSIZE=730,YSIZE=722



END
